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A century ago, the foundations of equilibrium statistical mechanics were laid. For a system in 

equilibrium with a thermal bath, much is understood through the Boltzmann factor, e , 

for the probability of finding the system in any microscopic configuration C. In contrast, apart 

*vj ■ from some special cases, little is known about the corresponding probabilities, if the same system 

is in contact with more than one reservoir of energy, so that, even in stationary states, there is a 
constant energy flux through our system. These non-equilibrium steady states display many surprising 

i-^ properties. In particular, even the simplest generalization of the Ising model offers a wealth of 

unexpected phenomena. Mostly discovered through Monte Carlo simulations, some of the novel 
properties are understood while many remain unexplained. A brief review and some recent results 
will be presented, highlighting the sharp contrasts between the equilibrium Ising system and this 
non-equilibrium counterpart. 
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I. INTRODUCTION 

>\ 

As we celebrate the centennial of the American Physical Society, we honor the founding of equilibrium statistical 
mechanics, which also took place about a century ago. That breakthrough enables us to understand thermodynamics 
}*J in terms of microscopies. Further, predictions based on the Boltzmann-Gibbs framework have been applied with so 
much success that we now take for granted many of the inventions of the industrial revolution, e.g., automobiles, 747's, 
power stations, etc. Yet, it may be argued that no systems are truly "in equilibrium", since infinite times and infinite 
^\ , thermal reservoirs or perfect insulations would be necessary. Indeed, essentially all natural phenomena bear the marks 
^^ ■ of non-equilibrium processes. Unlike the aforementioned class of "artificial" systems, most natural systems are not 
"set up" with special conditions, under which equilibrium statistical mechanics provides excellent approximations. 
H Unfortunately, the theory of non- equilibrium statistical mechanics is far less developed than its equilibrium counter- 

part. As a result, the most ubiquitous phenomena are the most poorly understood. In fact, relying on the intuition 
from equilibrium physics, we are often surprised, even by the behavior of systems in non- equilibrium steady states. 
These form a small subset of non-equlibrium phenomena where the states are time-independent, mimicking systems 
in equilibrium. In this article, the main differences between systems in equilibrium and non-equilibrium steady states 
will be highlighted. For example, stationary distributions of the former are well known. In contrast, we have no simple 
Boltzmann-like factor, e~ ' bT , f or n0 n-equilibrium steady states. Fortunately, with the aid of modern computers, 
it is possible to explore the behavior of simple model systems in stationary states far from equilibrium. An excellent 
example is the driven Ising lattice gas. Despite its simplicity, simulations continue to reveal a seemingly unending 
list of counter-intuitive phenomena. Yet, because of its simplicity, some of these surprises are now reasonably well 
understood. 

From textbooks, we learn that the first step in equilibrium statistical mechanics is to apply the fundamental 
hypothesis to a physical system in complete isolation: every configuration, C (or microstate), available to the system 
may be found with equal probability: i^ so [C] ex 1. By energy conservation, H[C] (the energy associated with C, 
which may include external static potentials) cannot change, i.e., A7Y =0 . Extending our scope to a system which 
can exchange energy with a much larger (theoretically infinte) reservoir and applying the fundamental hypothesis to 
the combination, we arrive at the canonical ensemble: when equilibriated, the probability for finding a system in C, 
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P eg [C], is proportional to exp(— H[C]/ksT), where T is the temperature associated with the thermal reservoir. In 
this stationary state, on the average, the energy of our system is a constant, which may be denoted by (ATi) T = 0. 
Fluctuations around this constant can be traced to losses or gains to the reservoir, while the average energy flux 
between them is zero. 

Here, we are interested in a system exchanging energies with two or more reservoirs, which are not coupled otherwise. 
If, say, the reservoirs are set at different temperatures initially, then we may expect the following scenario. Assuming 
the reservoirs are much larger than our system, then there should be a time when our system would be in a stationary 
state, while the reservoirs are still close to their initial states. In this sense, the combined system is far from equilibrium, 
with energy flowing from the hotter reservoir to the colder one. However, if we focus on our system alone, we find that 
its energy is constant on the average. Keeping in mind the presence of two reservoirs (for which we use the subscripts 
T and E), we denote this situation by (AH) T + (AH) E = 0. This state is also far from the equilibrium state, since 
there is a non-zero flux flowing through it. On the average, it gains a non-trivial amount of energy from one reservoir 
and loses energy to the other. In other words, neither term in the above equation vanishes: (AH) T = — (AH) E ^0. 
We will refer to such systems as being in non-equilibrium steady states. Examples abound in nature, from our planet 
as a whole to simple daily activities like cooking. The fundamental question for these states is: what is the stationary 
distribution P*[C]1 

Because our problem is inherently a time dependent one, we believe that the most appropriate approach is to start 
with the t-dependent distribution: P[C;t], the time evolution of which is governed by the master equation: 

^-P[C- i] = ]T {P (C -► C) P[C; t]-R(C^ C) P[C; t}} = LP. (1) 

m c 

Here R (C — > C) stands for the rate with which a configuration C changes to C and, in principle, can be found once 
we specify how our system is coupled to the various reservoirs. Then, the stationary P* will be "just" the (right) 
eigenvector of L with zero eigenvalue: = LP*. 

For physical systems, the task of finding the P's is clearly too complex. On the other hand, the success of 
equilibrium statistical mechanics suggests those P*'s are independent of the details of the rates. Apart from mixing 
(all the configurations being connected by the P's), the condition on the rates so that our system arrives at thermal 
equilibrium (being coupled to a single reservoir at temperature T) is known as detailed balance: 

X(C>^C)_(H[C]-H[m 



R(C->C) { k B T 

Numerous successful Monte Carlo simulations of systems in equilibrium are based on choosing the Metropolis rate 
0: R(C ->C) = min[l,e- (w[c ' 1 - w[cl)/fcsT ]. It is clear that, with rates satisfying eqn.(g), the Boltzmann P eq is a 
stationary distribution. Moreover, the equation = LP eg is satisfied by 

R (C -» C) P eq (C) = R(C^ C) P eq (C) 

for every pair (C ' ,C)\ 

An important distinction for a system evolving toward non-equilibrium steady states is that eqn.(0) no longer 
holds. The effects of two reservoirs cannot be embodied in eqn.(@). The stationary distribution P* is not known 
without solving = LP* first while each term on the right hand side of eqn.(nl) is not necessarily zero. A good 
analogue with electromagnetism is to regard eqn.(nj) as a continuity equation. Then each C is analogous to a node 
in a circuit while the terms on the right are net currents between pairs of nodes. While equilibrium corresponds to 
electrostatics (with time-independent charge distribution and zero currents), non-equilibrium steady states correspond 
to magnetostatics, for which currents are steady but non-vanishing. Given this sharp contrast between equilibrium 
and non-equilibrium steady states, it is not surprising that the latter problem is considerably more difficult, since 
P* itself is unknown a priori. On the other hand, the variety and richness is also much greater. In general, P* will 
depend on the details of the rates, although we expect the "universality classes" of P's leading to the same P* to 
be just as large as in the equilibrium case. In this short article, we will focus only on a particularly simple model - 
the driven Ising lattice gas (|]. Though the model appears simple, there is no analytic solution, so that its properties 
are generally explored through computer simulations. Now, most non-equilibrium systems are not translationally 
invariant. Typically, temperature gradients or velocity shears are present, so that the usual thermodynamic limits 
cannot be taken. As a result, the co-operative behavior in such systems is much more difficult to analyse than that in 
systems with translational invariance. An advantage of our simple model is that not only is it translationally invariant, 
it displays a host of surprising phenomena when driven to non-equilibrium steady states. In the next section, we will 
present a brief summary of the specifications of this model and some of the remarkable discoveries from Monte Carlo 



simulation studies. In many situations, the well honed arguments of equilibrium statistical mechanics, based on the 
competition between energy and entropy, fail dramatically. More details may be found in, e.g., J3[. Section 3 is 
devoted to some recent developments while some concluding remarks are included in the last. 

II. A BRIEF REVIEW OF THE DRIVEN LATTICE GAS 

Motivated by the physics of fast ionic conductors W, Katz, Lebowitz and Spohn introduced a simple model in 1983 
H, which has served as a primary testing ground for exploring unusual properties of non-equilibrium steady states. 
It consists of an Ising lattice gas |5|,|6| with attractive nearest-neighbor interactions, driven far from equilibrium by 
an external "electric" field, E. In the spin language, it is a ferromagnetic Ising model with biased spin-exchange M 
dynamics. Major reasons for choosing this model are: 

• having a translationally invariant dynamics, the steady state is expected to be also invariant; 

• the two reservoirs are coupled through an anisotropic dynamics; 

• many of its equilibrium properties are well-known, especially in two dimensions (d = 2) |8j,|9j ; 

• it is a system with non-trivial phases in d > 1, whether driven or not; 

• the equilibrium system can be reached continuously by taking the E — v limit; and 

• a different equilibrium system can be reached by choosing appropriate boundary conditions. 

For completeness, we give a brief description of the d = 2 model here. On a square lattice with fully periodic 
boundary conditions, each of the L x x L y sites may be occupied by a particle or left vacant, so that a configuration 
C of our system is completely specified by the occupation numbers {rii}, where i is a site label and n is either 1 or 
0. Translating the spin language of Ising Q is simple: s = In — 1 = ±1. An attractive interaction between pairs of 
particles in nearest-neighbor sites is modeled by the usual Hamiltonian: TL[C\ = — 4JX)<j ;> n i n ji with J > 0. The 
factor of 4 means that TL assumes the form — J^2, ss m the s pi n language, so that, in the thermodynamic limit, the 
system undergoes a second order phase transition at the Onsager g critical temperature Tq — (2.2692..) J/kg. For 
the lattice gas, this point can be reached only for half-filled systems, i.e., Y^i n i — L x L y /2. In Monte Carlo simulations 
for a lattice coupled to a thermal bath at temperature T, particles are allowed to hop to vacant nearest neighbor 
sites with probability min[l, e~ An ' kBT ] fil, where A7i is the change in Ti after the particle-hole exchange. Note that 
these rules conserve the total particle number ^V m, so that half-filled lattices must be used, if critical behavior is to 
be studied. Starting from some initial state, this dynamics should bring the system into the equilibrium state with 
stationary distribution P ce [C] oc e^ 7i ^ c ^ kBT . 

The deceptively simple modification introduced by Katz, et. al. 0] is an external "electric" field. Pretending 
the particles are "charged", the effect of the drive is to modify the hopping rates to biased ones. Specifically, the 
jump probabilities are now min[l, e -( AW - c - E )/ fc -B T ] ; where e = (—1,0,1), for a particle attempting to hop (against, 
orthogonal to, along) the drive and E is the strength of the field. Note that, locally, the effect of the external field 
is identical to that due to gravity. Indeed, had we imposed "brick wall" boundary conditions (particles reflected at 
the boundary, comparable to a floor or a ceiling), this system would eventually settle into an equilibrium state, like 
gas molecules in a room on earth. Of course, the reason behind this outcome is that gravitation is a static potential 
and can be incorporated into Ti[C]. The "price" paid is the loss of translational invariance due to the presence of 
boundaries, leading to inhomogeneous particle densities. 

Returning to our model, in which periodic boundary conditions arc imposed, we see that translational invariance 
is completely restored so that, in the final steady state, the particle density is homogeneous, for all T above some 
finite critical T c . At the same time, a particle current will be present. For gravity, such a situation exists only in art 
pQ|. In physics, this situation can be realised only with an electric field. Thanks to Faraday, a constant electric field 
circulating around the surface of a cylinder can be set up by applying a linearly increasing, axial magnetic field. If 
the particles are charged, they will experience the same force everywhere on the cylinder. As they also lose energy 
to the thermal bath, a steady state with finite current can be established. With this possibility in mind, we will use 
the term "electric" field to describe the external drive and imagine our particles to be "charged" . We should caution 
the reader that, unlike real charges, our particles are not endowed with Coulomb interactions between them, just like 
the typical neglect of gravitational attraction between gas molecules. Finally, note that a well defined, single- valued, 
potential which gives rise to such an electric field is necessarily time dependent. 

Given the microscopic model, we can ask: What are its collective properties when it settles down in a non-equilibrium 
steady state? In particular, what is P*[C]1 Since there is no global Hamiltonian, we cannot exploit Boltzmann's result. 
Instead, we must resort to the master equation (II]) and attempt to find the solution to = LP* . Though L is a sparse 
matrix, finding P* is a non-trivial task. Only for very small systems can P* be obtained [O. Of course, it is difficult 
to discern collective behavior such as phase transitions in "microscopic" systems like these. Nevertheless, we can 
already see that the P*'s here (fig. la) are quite distinct from P eq (fig. lb). Also, at this level, we can derive an 



interesting consequence: the violation of the standard fluctuation dissipation theorem. Computing 1(H) and (T~C 2 ), it 

is easy to verify that, in the driven case, (H 2 ) — (H) ^ —dp (Ti) . 

Turning to collective phenomena on the macroscopic scale, we face serious difficulties in finding P* analytically, 
let alone solving for thermodynamic quantities a la Onsager. Without modern computers, it would be impossible to 
make much progress. Using simulation techniques, we may answer questions like: what happens to the second order 
phase transition found by Onsager? In particular, how does the critical temperature depend on the drive, i.e., what 
is T C (E)1 Several simple possibilities come to mind: 

a) T C (E) jumps to infinity for any E , i.e., the drive, however small, orders the system. 

b) T C (E) rises with E indefinitely. 

c) T C (E) rises with E, saturating at some finite temperature T c (oo). 

d) T C (E) is independent of E , i.e., T C {E) = T . 

e) T C (E) decreases with E, either dropping to zero at finite E, or saturating at a finite T c (oo). 

f) T C {E) jumps to zero for any E , i.e., the drive, however small, disorders the system. 

While possibilities (a), (b), and (f) appear incredulous, intuitive arguments might be made for (d) and (e). The naive 
argument for (d) is that, in an inertial frame where the global current vanishes, the system should look just like an 
equilibrium Ising model. On the other hand, to arrive at (e), we think of the drive as gravity, feeding energy into the 
system, so that its effects should be the same as a reservoir with a temperature higher than the surrounding bath. 
Therefore, to order the system, the bath temperature would have to lowered, i.e., T C (E) < To- In reality, simulations 
P| offered the first surprise: T C (E) increasing with E and saturating at T c (oo) ~ lATo p2| ! A number of arguments 
now exists for this behavior, but all are "post-dictions" . Indeed, none of these are convincing, while approximate 
schemes for analytic computations of T C (E) offer only hints to the puzzle of T C (E) > T c (0) OJf4|. 

As the system is probed deeper, more surprises appear. Due to space limitations, we only list some of them here, 
refering the interested reader to PUB] for further details. 
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FIG. 1. Stationary distributions for a 4x2 Ising lattice gas, driven with infinite E (a) and in equilibrium (b). The normal- 
ization is set via: largest P* = 1. Corresponding configurations (only one out of the equivalence class) are shown, with the 
drive downwards. 



A. Disordered phase (T > T c ) 



In the equilibrium system, there is little of interest far above criticality. Correlations are short ranged so that 
most properties can be understood through Landau-Ginzburg mean field theory. When driven, however, this system 
displays 

• long range two-point correlations, decaying as \/r d [flql 

• singular structure factors, with a discontinuity at the origin @Jlj] 

• non-trivial three point correlations, the Fourier transforms of which show infinite discontinuities at the origin JL8[ ; 

• a fixed line, rather than the Gaussian fixed point, governing large-scale, long time behavior fll9[ ; 

• shape-dependent thermodynamics M 



B. Critical Behavior (T ~ T c ) 

In 1944, Onsager solved the d = 2 Ising model and computed many of its critical properties. However, the deeper 
understanding of critical phenomena came only in the 70's, with the advent of field theoretic renormalization group 
analysis pM. Within this framework, we learnt that a large class of systems fall into the Ising universality class, 
controlled by the Wilson-Fisher fixed point |22j. When driven into non-equilibrium states, 

only ordering into strips parallel to E occurs; 

only one of the two lowest structure factors diverges; 

strong anisotropy appears while longitudinal and transverse momenta scale differently; 

the critical dimension is 5 instead of 4; 

a new, non-Hamiltonian, fixed point and universality class is identified |23| ; 

a host of new exponents, though only one independent, e merg es; and 

anisotropic finite size scaling is essential for data collapse jl2| . 



We should note that, unlike properties far from T^, the critical properties were predicted by theory |23| well before 
confirmations from computer simulation studies |12| . 

C. Ordered phase (T <T C ) 

Below the critical point, phase segregation and co-existence occurs, with interfaces separating the particle rich from 
the particle poor regions. Due to the periodic boundary conditions, each region is a single strip wrapped around the 
torus. Unlike in equilibrium, only strips parallel to the drive exist. Correlations within each region are also expected 
to be long ranged, though they have not been studied carefully so far. Some intriguing properties of the interfaces, 
dramatically different from equilibrium 124] ones, are: 

• statistical widths remaining finite as L — > oo |25|| rather than diverging as vl; 

• the structure factor diverging as q~ - 67 |g6[ instead of the usual q~ 2 of capillary waves |27|]; • interface orientation 
affecting bulk energies pS] ; and 

• instabilities when forced to lie at a non- vanishing angle with respect to the drive |29(| . 
In addition to interface anomalies, other remarkable features include 

• coarsening during a quench showing several time regimes and asymmetries which cannot be accounted for by a 
modified Cahn-Hilliard theory [B0[ ; 

• systems subjected to shifted periodic boundary conditions displaying new, multistrip phases p8[ ; and 

• systems subjected to open periodic boundary conditions displaying icicle like fingers pi] . 

As we try to convey in this section, there are still many unresolved mysteries associated with this decepetively 
simple model. In the next section, we will present some recent, on-going investigations which attempt to probe deeper 
into the simple driven lattice gas. 

III. SOME RECENT DEVELOPMENTS 

A. Steady State Energy Fluxes 

In the introduction, we emphasized the importance of energy flow through our system in a non-equilibrium steady 
state. There has been no systematic study of this aspect, even though it may contain a key to the understanding of 
such systems. We begin by simply confirming our intuitive picture: that typical jumps parallel/trasverse to the field 
are associated with energy gain/loss. In the previous section, we showed the complete solution for a 4 x 2 lattice. 
Using those P*'s, we can compute these gains/losses and then compare them with Monte Carlo simulations |p2| . 
The results are best displayed as histograms for all possible values of AH after an attempted jump. There are two 
such histograms, associated with the two types of jumps. We also performed the same analysis for the equilibrium 
case. It is reassuring that simulations confirm, within statistical errors, all theoretical predictions. Encouraged by 
these results, we carried out simulations on a 30 x 30 lattice |32|] . Figure 2a shows that, for the equilibrium case, 
both histograms are entirely symmetric (within statistical errors), so that the gains/losses balance for cither types of 
jumps. By contrast, figure 2b shows asymmetric histograms, confirming our expectation that the systems tends to 
gain energy when a particle jumps in the field direction, etc. Our hope is that a good approximation scheme, within 
the theoretical framework described above, can be found leading to quantitative predictions of these histograms. 



B. A Possible New Phase in Lattice Gases with Anisotropic Interactions 

Since the drive introduces a non-trivial anisotropy into the Ising lattice gas, it is natural to ask how E might compete 
with anisotropic couplings. The equilibrium model is part of Onsager's solution, so that we may again compare driven 
cases with well-known results. Since the drive enhances longitudinal correlations, we had expected that U T C (E) will 
be higher (lower) if the drive is aligned with the stronger (weaker) bonds" |33|. Subsequent simulations showed the 
opposite |34|]! Indeed, with saturation drive, T C (E) can drop below T c (0) for a > 1.7, where a 2 is the ratio of the 
coupling along the field direction to the "transverse" coupling. Motivated to look into the behavior displayed at 
various temperatures, we find that the typical configurations are indeed disordered for T 3> T c (0) and ordered into 
a single strip for T <C T C (E). However, there is a significant range of temperatures where the system appears to be 
"ordered" in the drive direction without being in a single strip. In other words, while the densities in each column 
(i.e., along the drive) are bi-modally distributed pq] , the usual order parameter (5(1,0), structure factor with lowest 
transverse wave number) is still quite small. In more picturesque language, we call this a "stringy" state. Actually, 
this type of configurations have been previously reported |3q ] . But they were believed to be long transients on the 
way from disordered initial states to ordered ones and thus, disregarded. By contrast, we observed that, starting from 
ordered states, the system evolves towards, and spends considerable periods of time in, the stringy states. 
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FIG. 2. Histograms of ATI after an attempted jump in a 30 x 30 lattice, at T — 2.5 with E — (a) and E — oo (b). 

In an effort to quantify these states, we define the ratio: 

■ d G(0,L/2) 



n = L a 



5(1,0) 



(3) 



where G(x,y) is the (untruncated!) two point correlation function, with y being parallel to the drive, and 
S(k,p) = ^2 X G(x,y)e 27rilykx+pv ^ L is the structure factor, used ordinarily as the order parameter. In particular, 
we are interested in the behavior of 1Z as L — > oo . 

Far in the disordered phase, as a result of the r~ d decay, we have G(0,L/2) — > 0(L~ d ) while 5(1,0) — > 0(1). On 
the other hand, deep in the ordered phase G — > O(l) while 5 — * 0{L d ). Therefore, as long as the system is far from 
criticality, we have 



ft-vO(l). 



(4) 



Near criticality, if there are no stringy states, we may apply the results of renormalization group analysis p3| - 
G(0,r) -> r -[^-2+A]/(i+A) and s(k,0) -> k- 2 - to arrive at K -> z,M-3]A/(i+A)_ Here, A is the exponent associated 

with anisotropic scaling of the momenta (ku ~ fcj_ )■> an d is found to be (8 — rf)/3 to all orders in the expansion 

around (5 — d) > 0. Thus, for our simulations in d = 2, we see that 1Z decreases with L : 



K - L- 2 ' 3 



(5) 



By contrast, in a stringy state, ordering has set in for the longitudinal direction so that G(0, L/2) — > 0(1). Meanwhile, 
complete phase segregation into just two regions is yet to take place, so that 5(1,0) — > O(l). As a result, we expect 



K^L a 



(6) 



which is an increasing function of L. The differences between the three behaviors (Jj, |5|, and |3|) should be dramatic 
enough to discern in simulations. 

To provide further contrast, we perform the above analysis for Ising models in equilibrium. Above criticality, 
G — » e~ r /^r 2 ~ d and S — > 0(1) so that 1Z decreases exponentially as L 2 e~ L '^. Near T c , we have G — > r 2 ~ d+v and 
5 — > k~ 2+v leading to 7£ — > O(l). The behavior deep in the ordered phase is unchanged. So, in all cases, 1Z should 
not increase with L. 

Turning to simulations, we find that the stringy state seems to be most pronounced for large a. In particular, we 
focused on models with a = 3 and saturation drive and compiled data from systems sizes L = 10,20,30,40,60, and 
90, using two or three independent runs. Starting with both random and ordered initial conditions, the runs last up 
to 800,000 Monte Carlo steps. The plot of 1Z vs. L 2 in figure 3 shows that, for two of the temperatures investigated, 
this ratio appears to be increasing with L. Meanwhile, for T (fc [0.8, 1.0], the behavior is consistent with 1Z — > 0(1). 
These preliminary results lead us to conjecture the existence of a "stringy phase" , especially in the thermodynamic 
limit of square (L x L) lattices. Of course, simulations with larger L's will be needed to determine if the increasing 
behavior persists. Other systematic methods can also be brought to bear, such as distribution functions of both G and 
S. We have initiated a study of the histograms of S(k,p) for a range of low lying wave- vectors, using the time-series 
of each quantity. Preliminary data show remarkable structures which are being verified in longer runs. Attempts 
at a theoretical understanding, based on phenomenological approaches, of the nature of the "stringy phase" are in 
progress. 
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FIG. 3. Correlation- Structure Factor ratio 1Z for systems with a = 3 and saturation drive, plotted against L 2 . The legend 
refers to the temperatures of the various runs. 



IV. CONCLUDING REMARKS 



In this brief article, we highlighted several major differences between a system in thermal equilibrium and one 
in non-equilibrium steady states. Apart from the obvious presence of non-trivial energy fluxes, the latter systems 
display many distinguishing and surprising features. Since their stationary distributions are neither a priori known 
nor susceptible to analytic probes (except for some simple limits or 1-D cases), all efforts to uncover the macroscopic, 
collective properties of these systems are based on Monte Carlo simulations. 

Focusing on a particularly simple model - the Ising lattice gas, driven into non-equilibrium steady states by an 
external "electric" field 0, we gave a brief review of a variety of surprising and counter-intuitive behavior. In the last 
section, we offered two of the recent developments in this continuing saga: detailed investigations of the energy flux 



and preliminary studies indicating the possible existence of a new phase (especially for systems with large anisotropics 
in the attractive interactions) . Beyond this simple model, many generalizations have been explored. Examples include 
repulsive interactions, random drives, quenched random impurities, multilayers, and multispecies, to name but a few 
PL Further from this class of "driven diffusive systems" is a wide range of other non-equilibrium steady states, e.g., 
surface growth, electrophoresis and sedimentation, granular and traffic flow, biological and geological systems, etc. 
All of these offer further surprises, some understood and most unexplained. At present, each non-equilibrium system 
is studied independently from the others. The hope is that a unifying concept and framework, like the fundamental 
hypothesis or the Boltzmann factor, will be discovered before the next centennial meeting of the APS. 
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